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Abstract The Ching Hai Toad-headed Agama (Phrynocephalus vlangalii) complex is a small toad-headed viviparous 
lizard that is endemic to the Qinghai-Tibetan Plateau. A fragment of mtDNA ND4-tRNA""" from 189 samples in 
26 populations was used to infer the phylogeographic history of this species complex in the upper reaches of the 
Yellow River. Phylogenetic analyses revealed that P. vlangalii and another proposed species (Р. putjatia) do not 
form a monophyletic mtDNA clade, which in contrast with a previous study, includes P. theobaldi and P. forsythii. 
Lineage diversification occurred in the Middle Pleistocene for P. vlangali (ca. 0.95 Ma) and in the Early Pleistocene 
for P. putjatia (ca. 1.78 Ma). The uplift of the Anyemaqen Mountains and glaciations since the mid-late Pleistocene, 
especially during the Kunlun Glaciation, are considered to have promoted the allopartric divergence of P. vlangalii. 
The diversification of P. putjatia may be triggered by the tectonic movement in the Huangshui River valley during the 
C phase of Qingzang Movement. Subsequently, the glacial climate throughout the Pleistocene may have continued 
to impede the gene flow of P. putjatia, eventually resulting in the genetic divergence of P. putjatia in the allopatric 
regions. Demographic estimates revealed weak population expansion in one lineage of P. vlangalii (А2, the Qaidam 
Basin lineage) and one lineage of P. putjatia (B2, the north Qinghai Lake lineage) after approximately 42 000 years 
before present. However, constant population size through time was inferred for two lineages (Al and ВІ), the source 
of Yellow River lineage of P. vlangalii and the southeast of Qinghai Lake lineage of P. putjatia, possibly due to stable 
populations persisting in areas unaffected by glacial advances. Our results also suggest: 1) at least four differentiated 
lineages of P. vlangalii complex may have evolved allopatrically in different regions during the Pleistocene glaciation 
events; 2) in support of several recent studies, P. putjatia is a valid species, having a more wide distribution than 
previously considered; and 3) a hypothesis referring to P. у. hongyuanensis, inhabiting in the source region of the Yellow 
River, being synonymous with Р v. pylozwi is supported. 
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1. Introduction 


A fundamental premise of phylogeography is that 
geological events, climatic history, and environmental 
heterogeneity play an important role in population 
divergence and speciation (Avise, 2000). Phylogeographic 
analysis of DNA sequence data has been routinely used 
to determine species boundaries and speciation patterns, 
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especially in the context of Quaternary and Holocene 
histories of regional biotas (Hewitt, 2001). The dramatic 
geological and climatic history and striking habitat 
diversity of the Qinghai-Tibetan Plateau (QTP), ranging 
from broad plantation surfaces, basins to mountain ranges, 
have made it a geographic focus of many phylogeographic 
studies (Jin et al., 2008; Qu et al., 2010; Xu et al., 2010; 
Zhan et al., 2011). As reviewed by Yang et al. (2009), the 
geomorphic and climatic changes on the plateau during 
the Quaternary have a remarkable influence on regional 
and adjacent biogeographic patterns. 

The Quaternary glaciations in the QTP and the 
bordering mountains were the consequence of a 
combination between climate and local tectonic uplift. 
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In particular, the Kunlun-Huanghe and Gonghe tectonic 
uplifts have played very important roles in triggering 
glaciations in high Asia (Shi et al., 2011; Zhao et al., 
2011; Zhou et al., 2011). Although there is no evidence 
of a continuous ice sheet covering the whole Tibetan 
Plateau throughout the Quaternary glaciations (Li et al., 
1991; Zheng and Rutter, 1998; Zhou et al., 2011), four 
discrete Pleistocene glaciations have been identified on 
the QTP and the bordering mountains. These four main 
Pleistocene glaciations are correlated with MIS 18-16, 
12, 6 and 4—2 (Zhou et al., 2011; and references therein). 
The most extensive glaciation in the plateau, Kunlun 
Glaciation, occurred during MIS 18—16 (0.78—0.62 
Ma), which corresponds with the late Kunlun-Huanghe 
movement (0.7—1.1Ma, Cui et al., 1998). Several 
researchers suggested that current populations of some 
species originated from ancestors that recolonized this 
area from low-altitude refugia located at the plateau-edge 
during the Last Glacial Maximum (LGM, approximately 
20 000 years ago). However, if no massive ice sheet 
developed on the QTP during the Pleistocene glaciations 
(Zhou et al., 2011), it 1s highly likely that some cold- 
tolerant reptiles could have survived in ice-free areas of 
the plateau, even during glaciations. 

The Yellow River originates from the Yuegouzonglie 
Basin on the north foot of the Bayan Har Mountains 
of the QTP, winding eastward and crossing the Loess 
Plateau and the Huang-Huai-Hai Plain into the Bohai 
Sea. Its formation and evolution is closely related to the 
uplift of the QTP during the Quaternary period (Jiang 
et al., 2007; Li, 1993; Li et al., 1996; Liu and Sun, 
2007; Zhang et al., 2004). Li et al. (1996) discussed 
the geomorphological evolution of the Yellow River in 
the upper reaches and the relationship with QTP uplift, 
suggesting the initiation of the upper Yellow River 
was at 1.6 Ma near Lanzhou, and has experienced a 
strong incision since 0.15 Ma. Pan (1994) discussed the 
geomorphic evolution and development of the upper 
reaches of Yellow River in Guide Basin. He suggested 
that the Gonghe movement (0.15 Ma) made Guide Basin 
rise. The uplift and dry climate made Guide Basin lose 
the water link with the Qinghai Lake (also known as 
Koko Nor in Mongolian, and mTsho sngon po in Tibetan) 
in Late Pleistocene, and instigated the disappearance 
of the Gonghe-Guide-Qinghai palaeo-lake (Pan, 1994; 
Yuan et al., 1990). Uplifting mountain ranges crossing 
the river course would have acted as boundaries for the 
eastern Qinghai palaeo-lake, near the Longyang gorges 
(Gonghe Basin exit), Shanba gorges (Guide Basin exit) or 
Jishi Gorges (Xunhua Basin exit). Madsen et al. (2008) 


Xianguang GUO et al. Phylogeographic Analyses of the Phrynocephalus vlangalii Species Complex 53 


proposed that Qinghai Lake began to form during the 
Middle or Late Pleistocene when the Riyue Shan uplift 
separated the Qinghai Basin from the Yellow River 
drainage system. The occurrence and development of the 
Quaternary glaciers in the upper reaches of Yellow River 
were also closely related to the amplitude and time of 
mountain uplifting and climatic conditions (Deng et al., 
2004; Zheng and Wang, 1996). Since the mid-Pleistocene 
epoch, three glaciations have occurred in the Bayun Har, 
Zanjia, Buqingshan and A'nyemaqen Mountains in the 
upper reaches of the Yellow River. 

It has been long recognized that the phylogeny and 
taxonomy of the toad-headed agamas (Phrynocephalus) 
is very complex and controversial. Among the many 
controversies, the classification of the Ching Hai Toad- 
headed Agama (Phrynocephalus vlangalii) complex 1s 
perhaps one of the most confusing (Jin et al., 2008). 
Historically more than ten names have been applied to 
P. vlangalii (Barabanov and Ananjeva, 2007). As for P. 
putjatia, Wang et al. (2002) considered it a valid species 
based on morphological and chromosomal differences, 
whereas most herpetologists regard it as a synonym of Р 
vlangalii (Zhao, 1999). Guo and Wang (2007) suggested 
that P. putjatia 15 a valid species, comprising populations 
from Guide, Qinghai Province and Tianzhu, Gansu 
Province; while P. hongyuanensis is not a valid species, 
synonymized instead with P. vlangalii. Subsequently, 
Jin et al. (2008) implied that P. putjatia is valid, and P. 
v. hongyuanensis is the synonym of P. v. pylzowi on the 
basis of mtDNA ND2 fragment. Recently, Nobel et al. 
(2010) reconfirmed the validity of P. putjatia by using 
Bayesian model-based assignment tests of microsatellite 
DNA data. As for the distribution of P. putjatia, Bedriaga 
(1909) originally thought Guide County, Qinghai 
Province, China as the only known place. Until recently, 
P. putjatia was thought to occur in southeast of the 
Qinghai Lake (Jin, 2008). Accordingly, the distribution 
information for Р putjatia is still incomplete. 

Although all mtDNA genes belong to a single linkage 
group and evolve as one locus, different mitochondrial 
genes may have different phylogenetic performance 
(Mueller, 2006; Zardoya and Meyer, 1996). In addition, 
individual mitochondrial genes might potentially produce 
misleading evolutionary inferences and hence might not 
constitute an adequate representation neither of the entire 
mitochondrial genome nor of the evolutionary history of 
the organisms from which they are derived. Accordingly, 
it is necessary to test the hypothesis of Jin et al. (2008) 
to clarify the evolutionary history of the P. valangalii 
species complex. 
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Here, we examine the phylogeography of the P. 
vlangalii species complex from the upper reaches of the 
Yellow River by conducting molecular phylogenetic, 
dating and population genetic analyses of mtDNA 
ND4-tRNA'*", We also include several haplotypes 
of P. theobaldi and P. forsythii (Pang et al., 2003) as 
outgroup taxa to test the phylogenetic relationships 
among the viviparous group. This provides a framework 
for achieving the primary goal of testing relationships 
between physical events and/or climatic changes in the 
upper reaches of the Yellow River and genetic diversity 
among populations. The major goals of this study 
were to address the following questions: 1) What is 
the geographic structure of morphological and genetic 
variation among populations of the upper reaches of the 
Yellow River? 2) Did the uplifting of the A'nyemagen 
Mountains or past glaciations promote the allopartric 
divergence of P. vlangalii and if so, when did it happen? 
3) Was there a single refugium or multiple refugia in the 
upper reaches of the Yellow River and adjacent areas? 
In addition, we discuss the taxonomic implications of 
our preferred hypothesis relative to previous notions of 
viviparous Phrynocephalus species relationships (Guo 
and Wang, 2007; Jin et al., 2008; Pang et al., 2003). 
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2. Materials and Methods 


2.1 Geographic sampling One hundred and seventy five 
Phrynocephalus vlangalii species complex specimens 
were collected from 24 sites across nearly the upper 
reaches of the Yellow River (Figure 1; Table 1). After 
being euthanized, liver or muscle tissues of individuals 
were dissected and preserved in 95% ethanol for 
genomic DNA extraction. Voucher specimens are held 
at the Chengdu Institute of Biology, Chinese Academy 
of Sciences. 14 published sequences (Table 1; Pang et al., 
2003) were also included in this study. Thus, a total 
of 189 individuals from 26 sites were included in this 
study. Voucher specimens and geographic localities 
corresponding to tissues used in this study are detailed in 
Table 1. Р axillaris and four other species were selected 
as outgroup taxa based on current understanding of the 
phylogenetic relationships among toad-headed lizards 
(Guo and Wang, 2007). The 16 sequences of the outgroup 
taxa were all retrieved from the GenBank. 


2.2 DNA extraction, amplification and sequencing 
protocols Total genomic DNA was extracted from 
fixed liver or muscle tissue following the method 
of Aljanabi and Martinez (1997). A fragment of the 
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Figure 1 A map of the upper reaches of the Yellow River showing the sampling sites of P. vlangalii species complex (numbered red dots). 
Site names are listed in Table 1. Distributions of the main lineages identified in the phylogenetic trees are shown. A: Riyue Shan; Ж: Xining 


City; [0: Huzhu County. 
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mitochondrial gene encoding the fourth subunit of NADH 
dehydrogenase (plus downstream Serine, Histidine, and 
Leucine tRNAs; hereafter collectively referred to as 
ND4-tRNA“”) was amplified via PCR using the primers 
described in Arévalo et al. (1994). Amplification was 
performed on the PTC-200 thermal cycler (Bio Rad) in 
a total volume of 50 uL containing 10 х reaction buffer 
(10 mM Tris-HCl, pH 8.3, 50 mM KCI, 1.5 mM MgCl), 
0.1mM of each dNTP, 0.4 uM of each primer, 1.25 units 
rTaq polymerase (TaKaRa) and approximately 100 ng 
of template DNA. Thermal cycling was performed with 
an initial denaturation for 4 min at 94 °C followed by 35 
cycles of 40 s at 94 °C, 40 s at 55 °C, 1 min at 72 °C with 
a final extension of 7 min at 72 °C. Negative controls 
were run for all amplifications. The amplified products 
were purified on a 0.8% agarose gel stained with ethidium 
bromide, using the DNA gel extraction kit (Omega Bio- 
Tek). The purified products were sequenced directly 
by Invitrogen Trading (Shanghai) Co., Ltd., using the 
same primers as PCR. The sequences were deposited in 
GenBank with accession numbers EU294028-EU294095, 
JF736517—]F736619. 


2.3 Sequence alignment and analyses A total of 
205 sequences were first aligned using ClustalX 1.83 
(Thompson et al., 1997) with default gap penalties. These 
alignments were rechecked based on inferred amino 
acid sequence homology (protein coding region) and 
secondary structure models (for tRNA regions; Macey 
and Verma, 1997) using MEGA v.4 (Tamura et al., 2007). 
This final alignment consisted of a total of 703 aligned 
positions (bases), three of which were inferred as gaps 
(relative to outgroup sequences) occurring in the loop 
structures of tRNAs in P. vlangalii. The data matrices are 
available from the first author. 

Compositional heterogeneity was evaluated using 
Chi-square (y^) tests implemented in PAUP* 4.0b10 
(Swofford, 2003) and assessed using the software SeqVis 
v.1.5 (Ho et al., 2006) to visualize base composition 
and to conduct matched-pairs tests of symmetry of base 
substitution (Ababneh et al., 2006). Evidence of evolution 
under conditions more complex than that assumed by 
commonly applied models (i. e., stationary, reversible 
and homogeneous conditions) was inferred if the scatter 
of dots in the tetrahedral plots was widely dispersed and 
If a proportion X of the matched-pairs tests of symmetry 
rejected symmetry with P-values less than or equal to 
X. This procedure is consistent with that advocated by 
Jermiin et al. (2008, 2009). If the matched-pairs tests 
yield a larger than expected proportion of probabilities < 
0.05 (1. e., with > 596 of the tests producing probabilities 


Xianguang GUO et al. Phylogeographic Analyses of the PArynocephalus vlangalii Species Complex 59 


< 0.05), then the conclusion is that the sequences have not 
evolved under stationary, reversible and homogeneous 
conditions. 


2.4 Phylogenetic analyses Phylogenetic hypotheses 
were generated with ND4-tRNA“” segments using two 
commonly applied phylogenetic methods: heuristic 
searches using equally weighted maximum parsimony 
(MP) analyses performed with the program PAUP* and 
Bayesian inference (BI) with the program MrBayes v.3.2 
(Ronquist and Huelsenbeck, 2003). In both MP and BI 
analyses, each haplotype was treated as a taxon and each 
nucleotide was treated as a character, and gaps were 
treated as missing data. 

For heuristic searches under parsimony, invariant 
characters were removed from the dataset, and all 
remaining characters were treated as equally weighted. 
Each search involved ten random addition replicates, one 
tree held at each step, TBR branch swapping, steepest 
descent on, and a maximum of 100 000 saved trees; 
all other search settings were left at default values. 
Non-parametric bootstrapping was used to generate 
phylogeny confidence values (Felsenstein, 1985), with 
1 000 pseudoreplicates using a heuristic tree search for 
each pseudoreplicate. These are presented as a majority 
consensus tree where branches with support lower than 
50% are collapsed. P. axillaris was used to root the trees 
based on a recent publication (Guo and Wang, 2007). 
Because intraspecific gene evolution cannot always be 
represented by a bifurcating tree, haplotype networks 
may more effectively portray the relationships among 
haplotypes within species [reviewed by Posada and 
Crandal (2001)]. Therefore, we constructed unrooted 
parsimony networks of haplotypes for P. vlangalii species 
complex using TCS v.1.21 (Clement et al., 2000). 

The appropriate models for the BI analyses were 
selected using Akaike Information Criterion (AIC; 
Akaike, 1974) in the program MrModeltest 2.3 (Nylander, 
2008) with each starting tree obtained using the neighbor- 
joining algorithm in PAUP*. Sites in the ND4-tRNA"" 
alignment were partitioned a priori into classes based on 
general evolutionary constraints, and Bayes factors were 
calculated to determine the optimal partitioning scheme 
(Brandley et al., 2005). Two partitioning schemes were 
assessed: one partition (one model for all sites), four 
partitions (1* codons, 2" codons, 3" codons, all other 
tRNA sites). The marginal likelihood of each partitioning 
scheme was approximated using the stepping-stone 
approach (Fan et al., 2011; Xie et al., 2011) implemented 
in Phycas (Lewis et al., 2008). Bayes factors strongly 
favored four partitions over one partitioning scheme (2In 
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Bayes factor = 242.282). Models chosen were as follows: 
GTR + G + I (1* codon position); HKY (2™ codon 
position); GTR (3" codon position); HKY + I (tRNAs 
positions). 

We estimated posterior probability distributions by 
allowing four incrementally heated Markov chains (default 
heating values) to proceed for 2 х 10’ generations, with 
samples taken every | 000 generations. Analyses were 
repeated beginning with different starting trees to ensure 
that our analyses were not restricted from the global 
optimum (Huelsenbeck ег а!., 2002). MCMC convergence 
was explored by examining the potential scale reduction 
factor (PSRF; Gelman and Rubin, 1992) convergence 
diagnostics for all parameters in the model (provided by 
the sump and sumt commands) and graphically using 
the program Tracer v.1.5 (Rambaut and Drummond, 
2009). The first four million generations, before this 
chain reached apparent stationarity, were discarded, and 
the remaining samples from the independent runs were 
pooled to obtain the final approximation of the posterior 
distribution of trees. To yield a single hypothesis of 
phylogeny, the posterior distribution was summarized as a 
50% majority-rule consensus. 


2.5 Bayesian hypothesis testing Bayes factors were 
used to compare the preferred Bayesian tree topology (see 
below) to Bayesian trees with constraint. This method 
differs from traditional hypothesis testing because it 
does not offer a criterion for absolute rejection of a null 
hypothesis but instead an evaluation of the evidence in 
favor of the null hypothesis (Kass and Raftery, 1995). The 
phylogeny inferred from the ND4-tRNA“” data set was 
constrained to alternative hypotheses. Constraint analyses 
were conducted in MrBayes v.3.2 using the command 
prset topologypr = constraint. All analyses consisted 
of two simultaneous runs each with an abbreviated 
three MCMC chains run for 20 million generations. 
The Bayes factor was determined by calculating the 
marginal likelihood for both unconstrained and constraint 
analyses using TRACER v.1.5 (Rambaut and Drummond, 
2009). The difference in these In-transformed marginal 
likelihoods was compared to the table provided by 
Raftery (1996). Based on the table, we consider а 2In 
Bayes factor > 10 as significant evidence for a hypothesis 
(Kass and Raftery, 1995). 


2.6 Divergence-time estimation BEAST v.1.6.1 
(Drummond and Rambaut, 2007) was used to estimate 
divergence time of lineages under a Bayesian statistical 
framework. As no direct fossil records of Phrynocephalus 
have so far been detected in Chinese fossil layers 


(Zerova and Chkhikvadze, 1984), we used the average 
nucleotide substitution rate of ND4, 1.17% (with 95% 
credibility interval 0.85%—1.81%) per site per million 
years (estimated from the data of Guo and Wang, 2007) 
among lineages, to date lineage divergences based on the 
uncorrelated lognormal relaxed clock model and a Yule 
process tree prior (Drummond et al., 2006). Analyses 
were performed using an HKY+G substitution model 
(best-fit model based on AIC in MrModeltest). The 
MCMC chain was run for 2 x 10’ generations, and it 
was sampled every 1 000 generations with the first 10% 
discarded as burn-in. Convergence of the MCMC chain 
to stationary distribution and adequate sampling were 
checked in TRACER. We performed the analysis three 
times to confirm the stability and convergence of the 
MCMC chains. The results were found to be satisfactory 
and samples from the multiple runs were combined. 


2.7 Demographic analyses Genetic diversity estimates 
were calculated for separate lineage of the P. vlangalii 
species complex using DnaSP v.5.0 (Librado and 
Rozas, 2009). These include haplotype diversity (Л) and 
nucleotide diversity (л). Mismatch distributions were 
calculated for each diagnosed population lineage to 
infer changes in population size with ARLEQUIN v.3.5 
(Excoffier and Lischer, 2010). To compare the observed 
data with the expected data under the sudden expansion 
model, we conducted goodness-of-fit tests based on 
the sum of squared deviations (SSD) and Harpending's 
raggedness index (r,) (Harpending, 1994) using 10 000 
parametric bootstrap replicates. Although the raggedness 
index (r,) is often used to assess significance of mismatch 
plots, multi-modal distributions that fit sudden-expansion 
models indicate structuring within populations or 
populations that are stable or shrinking (Rogers and 
Harpening, 1992). Therefore we do not consider the good 
fit of sudden-expansion models to unimodal mismatch 
distributions strong evidence of population expansion. 
Given that tests based on the mismatch distribution 
have little power against population growth, Tajima's 
D (1989), Fu’s Е; (Fu, 1997) and R2 statistic (Ramos- 
Onsins and Rozas, 2002) were calculated in DnaSP as 
additional assessments of possible population expansion. 
Significance of D, F and R2 values was determined 
in DnaSP using 1 000 simulated samples to produce 
an expected distribution under selective neutrality and 
population equilibrium. The cut-off level for statistical 
significance was 0.05. For Fu's F, significance at the 0.05 
level was indicated when P values were « 0.02 (Excoffier 
and Lischer, 2010). If significant statistic-values of 
Fu's Fs, Tajima's D and R2 are obtained, a population 
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Table 1 The sampling locations, number of specimens (n), and voucher number for the samples of P. vlangalii complex used in this study. 


Site label Location code Sampling location n Coordinates Voucher numbers" 

1 НУ Hongyuan 9 N33.17649 E102.62830 DLW200607001—2, DL050002-5, 2LDW001, 0714, 0541 

2 XMa Xiaman A 17 N33.70605 E102.48526 DLW200607007—9, DLW200607011-15, 
DLW200607017—23, Xm276, Xm277 

3 XMb Хіатап B 9 N33.77507 E102.54921 DLW200607024—26, DLW200607029-34 

4 XMc Xiaman C 13 N33.72512 E102.46995 DLW200607035—40, DLW200607042-46, 
DLW200607048, DLW200607050 

5 XMd Xiaman D 8 N33.74796 E102.50430 DL050011—15, DL050017-18, DL050020 

6 MQa Maqu А 8 N33.95074 E102.08892 DLW200607052, DLW200607054—55, DLW200607057— 
58, DLW200607060—62 

7 МОБ Maqu В 8 N33.89865 E102.12673 DLW200607064—70, DLW200607072 

8 GDA Guide A 11 N36.08161 E101.46884 DLW200607074-76,78-83; 0771,0772 

9 GDB Guide B 8 N36.27295 E101.08978 DLW200607084-87,103—106 

10 GDC Guide C 7 N35.91614 E101.32405 2WDL0508028-30, 32, 33, 35, 38 

11 GNA Guinan A 5 N35.79170 Е101.13223 2WDL0508044—48 

12 GNB Guinan B 2 N35.50655 E101.10550 2WDL0508059, 60 

13 HM Heimahe 11 N36.65030 E100.01867 DLW200607138-148 

14 GHA Gonghe A 5 N36.15008 E100.92765 DLW200607109-113 

15 GHB Gonghe B 7 N36.19800 E100.74776 DLW200607114—120 

16 TZ Tianzhu 8 N37.11158 E102.99005 WGL070602—08; 0709, 0769 

17 HYA Haiyan A 8 N36.83891 E100.84207 DLW200607123-128,135 

18 HYB Haiyan B 6 N37.18976 E100.55337 DLW200607155—160 

19 GC Gangcha 6 N37.31397 E100.15899 DLW200607 149-154 

20 XHA Xinghai A 4 N36.01592 E100.25779 2WDL0508052, 53, 55, 56 

21 XHB Xinghai B 6 N35.91789 E100.03450 GXG07080071,74-78 

22 MD Maduo 8 N35.20496 E98.97128 GP0708073-80 

23 DL Dulan 6 N36.01934 E97.95179 0792, 0793, 0796, 0797; GP0708159, 162 

24 GM Golmud 6 N36.38700 E94.96739 GP0708255, 257, 258; Кад1–3 

25 AKS Akesai 1 N39.41941 E94.28142 740 

26 SGH Suganhu 2 N38.89919 E93.89908 0970, 0972 


': The DNA data of the samples in italic were retrieved from GenBank. 


expansion model cannot be rejected. 

To further evaluate the signature of population 
fluctuations through time, we constructed Bayesian 
skyline plots (BSP; Drummond et al. 2005) implemented 
in BEAST. As above, the appropriate model of nucleotide 
substitution for each lineage was determined using 
MrModeltest. Separate analyses were performed for each 
lineage, with the exception of A3. Multiple analyses were 
performed to test for convergence, with each data set 
run twice for 2 x 10’ iterations and sampled every 1 000 
iterations with the first 10% initial samples discarded as 
burn-in under a strict clock, stepwise skyline model and 
uniform priors. Results from replicate runs were pooled 
using LogCombiner and skyline plots generated using 
TRACER. To convert estimates into units of time, we 
used mean substitution rate of 1.596 per site per million 
years as a heuristic, which was obtained from the BEAST 
analysis (See above). 


3. Results 


3.1 Base composition and nucleotide substitution 


patterns The sequence comprised 703 aligned sites, 
including 578 bp of partial ND4, 125 bp of tRNAs for 
all samples of the P. vlangalii complex. Alignment 
was unambiguous and 3 indels were inferred among Р 
vlangalii species complex. No premature stop codons 
were observed in ND4, indicating that the obtained 
sequences were mitochondrial in origin and not nuclear 
pseudocopies. A total of 39 haplotypes were revealed 
within the 189 individuals by 110 polymorphic sites 
containing 4 singleton variable sites. No haplotypes were 
geographically widespread. Haplotypes shared across 
populations were just found in some geographically 
proximate localities and the rest 32 ones were private for 
each locality (Table 2). The majority of variables of the 
703 bp were from ND4 (85.5%), with the third codon 
position (53.6%). The first (24.5%) and second (7.3%) 
variable positions were relatively low. The average 
proportions of base composition were T = 0.252, C = 
0.261, А = 0.361, and С = 0.126. Nucleotide composition 
showed an anti-G bias and 4+7-гісһпеѕѕ (61.3%), which 
is a characteristic of the mitochondrial genome. 

A base stationarity test shows insignificant differences 
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Мо. 1 


among all taxa (including 39 haplotypes and 16 outgroup 
taxa sequences) in base composition bias in the ND4 
data (X = 55.64, df = 162, P = 1.00). However, base 
composition at the third position exhibited significant 
heterogeneity: first position, Х = 15.99, df = 162, P = 
1.00; second position, x^ = 2.97, df= 162, P = 1.00; third 
position, % = 209.26, df = 162, P = 0.007; and tRNAs, 
х = 8.89, df = 162, P = 1.00. When all the sites are 
considered equal (1. e., all the sites placed in the same 
bin) and the tetrahedron is allowed to rotate, the 55 points 
are scattered tightly in an area where the proportion of А 
and C are > 25% (Appendix 1A). The points are clearly 
spread within a confined area, implying that there may 
be more compositional heterogeneity in these data than 
the initial analysis suggested. To address this issue, 
we binned the nucleotides according to the structural 
information. The distributions of points differ for the 
first, second, and third codon sites, with first codon 
sites displaying a small amount of scatter (Appendix 
1B), second codon sites displaying hardly any scatter 
(Appendix 1C), and third codon sites displaying a lot of 
scatter (Appendix 1D). The tRNAs sites also displayed 
a small amount of scatter (Appendix 1E). Rotating the 
four tetrahedral plots shows that the centroids differ for 
the codon sites and tRNAs sites, thus suggesting that it 
would be necessary to apply four Markov models to these 
data to analyze them appropriately within a phylogenetic 
context. To corroborate whether this is the case, the 
matched-pairs test of symmetry is used in conjunction 
with ND4-tRNA'""", 17, 2" 3" codon sites and tRNAs sites 
of the alignment. Table 3 summarizes the distribution of 
P-values for ND4-tRNA""", 1“, 2" and 3" codon sites and 
tRNAs sites. Of the 1.081 pairwise tests conducted: 1) 
seventeen tests (0.0157) for 1* codon site were found to 
produce a probability < 0.05, implying first codon site is 
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consistent with evolution under stationary, reversible, and 
homogeneous conditions; 2) seven tests (0.0065) for 2" 
codon site were found to produce a probability < 0.05; 3) 
two hundred and thirty five (235) tests (0.2714) for third 
codon site were found to produce a probability < 0.05, 
implying 3" codon site is inconsistent with evolution 
under stationary, reversible, and homogeneous conditions. 
Given these results, a sensible approach to analyze these 
data phylogenetically would be to apply a time-reversible 
Markov model to the 1* codon sites, another such model 
to the 2" codon sites, and a more general Markov model 
to the 3" codon sites. 


3.2 Phylogenetic analyses The heuristic search of 
the ND4-tRNA""" matrix resulted іп 228 871 equally 
parsimonious trees of 423 steps, with moderate CI 
(0.5816) and high RI (0.9052). These are presented 
as a majority-rule consensus tree where branches 
with bootstrap support lower than 50% are collapsed 
(Appendix 2). For the BI analyses, the 50% majority 
consensus tree is shown in Figure 2. The average PSRF 
was 1.000, implying convergence of runs. The MP and 
Bayesian trees differed in bootstrap support values, 
but not in the compositions of the major lineages. The 
parsimony bootstrap support was also marked on the 
branches that receive such support in the Bayesian tree 
(Figure 2). The two analyses strongly support monophyly 
of the viviparous group with high support values (BP 
= 10095; PP = 1.0). Haplotypes corresponding to P. 
theobaldi and P. forsythia were found to be nested within 
the main Phrynocephalus vlangalii complex clade. 
Clade B (Р. putjatia) is the sister taxon to the remaining 
viviparous species. Accordingly, P. putjatia and the other 
putative subspecies of P. vlangalii together did not form 
a single monophyletic clade. These comprise two distinct 
clades (A and B), being further subdivided into smaller 


Table 3 Matched-pairs tests for stationarity of base composition. Using 1 081 tests for symmetry of base substitution, the number and 
proportion of tests rejecting the hypothesis of symmetry at the threshold level in column 1 are listed separately for tests using all ND4- 
tRNA‘ base sites, 17, 2" and 3" positions as well as tRNAs sites only. Because the proportions of statistically asymmetrical result of 3" 
codon positions exceeds the threshold levels, we judge the 3" codon positions primarily inconsistent with a model of base compositional 


stationarity. 
Threshold* ND4-tRNA'*" 1“ codon positions 2" codon positions 3" codon positions tRNAs 
Number Proportion Number Proportion Number Proportion Number Proportion Number Proportion 

0.05 229 0.2118 17 0.0157 7 0.0065 235 0.2714 0 0 
0.01 54 0.0499 0 0 0 0 65 0.0601 0 0 
0.005 32 0.0296 0 0 0 35 0.0324 0 0 
0.001 5 0.0046 0 0 0 0 4 0.0037 0 0 
0.0005 1 0.0009 0 0 0 0 0 0 0 0 
0.0001 1 0.0009 0 0 0 0 0 0 0 0 
0.00005 1 0.0009 0 0 0 0 0 0 0 0 


^: The smallest P value for ND4-tRNA""" is 0.0000, for the 1“ codon positions is 0.0455, for the 2" codon positions is 0.0237, for the 3" 


codon positions is 0.0005, and for tRNAs positions is 0.0719. 
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Figure 2 ND4-tRNA""" majority-rule consensus tree inferred from partitioned Bayesian inference by using MrBayes v.3.2, associations with 
less than 0.5 posterior probability were collapsed. Bayesian posterior probabilities, maximum parsimony bootstrap values are shown. 


subclades: (i) source of Yellow River populations (A1), 
Qaidam Basin populations (A2), and Suganhu population 
(A3); and (ii) populations inhabiting the Qinghai Lake 
basin, corresponding to P. putjatia (B), which are 
subdivided into populations found to the north (B1) and 
southeast (B2) of the lake. However, the relationships 
between A1, A2 and A3 are ambiguous. 

To get additional insight into the relationships among 
the P. vlangalii complex, we analyzed our data set, 
using the coalescent-based statistical parsimony network 
approach. The five unconnected networks revealed by TCS 
are given in Figure 3. These generally correspond to the 
five described phylogenetic lineages. The two networks 
corresponding to lineages Al and A2 are connected 


by H36 and H32, with eight mutation steps. There are 
several unobserved intermediate haplotypes, representing 
unsampled or extinct haplotypes. The haplotype network 
of lineage B1 showed a star-like clustering. H2 (from 
GuideA, site 8) seemed to be a central haplotype, with six 
other haplotypes connected by 1–3 mutation steps. 


3.3 Bayesian hypothesis testing Bayes factor analysis 
reflects our primary interests of evaluating the alternative 
hypotheses of interrelationships of viviparous species 
of Phrynocephalus. As mentioned above, analyses of 
the ND4-tRNA""" data yielded a non-monophyletic P. 
vlangalii complex. Bayes factor analyses of the ND4- 
ТЕМА" data were conducted to compare topologies with 
constraints to the Bayesian tree topology. There was very 


Мо. 1 


(A) Clade А 
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Source Region of Yellow River 


О © 
H39 H35 
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Southeast of 
Qinghai Lake 
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Figure 3 Statistical parsimony network showing genetic 
relationships and distance among 39 haplotypes of P. vlangalii 
species complex from different populations. In the network, hollow 
circles indicate sampled haplotypes, and small solid circles indicate 
unsampled or extinct haplotypes. Each mutation step is shown as 
either a short or longer line connecting neighboring haplotypes 
(including observed and unobserved one). The size of the hollow 
circles roughly represents the numbers of samples carrying the 
haplotype, with the scale given beside the network. 
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strong evidence against the constrained monophyly of P 
vlangalii complex (Ln Lynconstrained = —3033.657, [01 constrained 
= —3067.845; 211 Bayes factor = 28.376, > 10). Thus, 
the monophyly of P. vlangalii complex is significantly 
rejected. 


3.4 Divergence times Bayesian estimates of divergence 
times are given for the major nodes in Figure 4. Dates 
ranged from 4.84 Ma (95% highest posterior density 
interval-HPD of 2.97—8.03 Ma) for the initial divergence 
of viviparous phryncoepahlus group to 0.95 Ma (95% 
HPD: 0.48–1.71 Ма) for cladogenesis within clade A and 
0.29 Ma (95% HPD: 0.09–0.65 Ma) for divergence within 
lineage B2. 


3.5 Historical demography Haplotype diversity (/) and 
nucleotide diversity (z) within each lineage were listed 
in Table 4. Genetic diversity of lineage B2 was lower 
than those of other lineages. Among four lineages, the 
highest haplotype diversity and nucleotide diversity was 
from lineage ВІ (A = 0.892, л = 0.0087). Graphs of the 
mismatch distribution (Figure 5) indicate a significant 
unimodal distribution suggesting expansive population 
growth for all lineages, although the lineage A2 appears 
slightly bimodal. Confirming the mismatch distribution, 
tests of the observed SSD and 7, statistics in these lineages 
were small and not statistically significant at the 0.05 level 
(Table 4), indicating that the sudden expansion model 
could not be rejected for any of the lineages. However, 
this pattern was not detected by Tajima's D, Fu’s Fs and 
R2 statistics, with non-significant statistic values for all 
lineages (Table 4). This further indicated that all lineages 
were evolving in a neutral manner, i.e., in mutation-drift 
equilibrium. This incongruence can be explained by the 
fact that tests based on the mismatch distribution have 
little power against population growth (Ramos-Onsins and 
Rozas, 2002). It is noted that the number of segregating 
sites has a substantial effect on the power of Tajima's 
D, Fu's Fs and R2 tests, which has been shown to have 


Table 4 Summary of genetic diversity indices, mismatch distribution analysis and neutrality tests. Number of individuals (Л), haplotypes (К) 
and segregating sites (S) , haplotype diversity (/), nucleotide diversity (л), sum of squared deviations (SSD), Harpending's raggedness index 


(ra), Tajima’s D, Fu's Fs and R2 tests are indicated. 


Lineage N K S h л та SSD R2 statistic  Tajima's D Fu's Fs 
AI 80 10 9 0.8351 0.0026 0.03469 0.00153 0.1013 0.00855 -1.513 

A2 23 6 8 0.7075 0.00293 0.16532 0.07852 0.1224 -0.1679 0.7109 
Bl 64 15 27 0.8924 0.00872 0.03623" 0.01936 0.1128 0.228 0.6273 
B2 20 6 5 0.6737 0.00155 0.04654 0.00037 0.1076 -0.7059 -2.016 

A1+A2+A3 105 17 31 0.8908 0.0083 0.01483 0.02076 0.0931 -0.05095 0.6983 
B1+B2 84 21 49 0.92 0.01934 0.0303 0.04242 0.1379 1.2424 3.2063 


The results of the observed mismatch distribution against a sudden expansion distribution are reported as the Harpending (1994) raggedness 
r, statistic and the sum of squared deviations (SSD). No mismatch distribution statistics were significantly different from the sudden 
expansion model, with exception to r, for lineage ВІ. However, no neutrality tests statistics were significant. АП values for each of these tests 


are non-significant above the threshold P = 0.05 except `. 
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Figure 4 Chronogram showing results from the uncorrelated lognormal Bayesian relaxed clock in BEAST v.1.6.1. Numbers and the numbers 
in brackets below the nodes are the posterior mean ages and 95% HPD, respectively. 


more power against population growth than tests based 
on the mismatch distribution such as r,. The number of 
segregating sites (S) ranged between 5 and 14 for the four 
lineages tested, which is below the point at which these 
three tests lose power. Low number of segregating sites 
also results in similar mismatch distributions for both 
constant populations and expanding populations (Ramos- 
Onsins and Rozas, 2002). For this reason, coalescent 
methods such as BSP, which take genealogy into account, 
may provide a better estimate of demographic history 


given these data. 

The coalescence-based Bayesian demographic analyses 
depicted the detailed fluctuations of effective population 
size for the four lineages (Figure 6). Median estimates 
from the BSP revealed weak population expansion in 
one lineage of P. vlangalii (lineage A2, the Qaidam 
Basin lineage) and one lineage of P. putjatia (lineage 
B2, the north of Qinghai Lake lineage). It should be 
noted that error around these estimates of Ne through 
time must be considered and this reduces certainty in 
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Figure 5 Mismatch distribution for observed frequencies of each lineage of P. vlangalii species complex compared to the simulated 


frequencies for sudden (stepwise) population expansion. 


demographic trajectories for all lineages. Lineages A2 and 
B2 appear to have experienced a long period of relatively 
constant population size, followed by weak expansion 
after approximately 42 000 years before present. The 
population growth rate of A2 was higher than that of B2. 
Compared with A2 and B2, lineages A1 (the source of 
Yellow River lineage) and ВІ (the southeast of Qinghai 
Lake lineage) appear to have remained stable over the last 
200 000 and 650 000 years, respectively. 


4. Discussion 
4.1 Phylogenetic relationships of Phrynocephalus 


vlangalii species complex Monophyly of the viviparous 
group was recovered from ND4-tRNA'""" data, which is 


consistent with the previous hypothesis (Guo and Wang, 
2007; Pang et al., 2003). To date, however, the intra- 
relationship of viviparous species is still equivocal. Guo 
and Wang (2007) proposed that the P. vlangalii complex 
is paraphyletic with P. theobaldi /P. zetangensis on the 
basis of analyzing 2065 bp mtDNA data (128, 168, cyt 
b, ND4-tRNA‘“" у, Later, Jin et al. (2008), on the basis 
of 609 bp mtDNA ND2-tRNA ^" fragments, revealed 
that P. vlangalii and two other species (P. erythrurus 
and P. putjatia) together form a clade which excludes 
P. theobaldi and P. zetangensis. The main discrepancy 
of the two hypotheses lies in the different phylogenetic 
placement of P. putjatia. In support of Guo and Wang 
(2007), our results corroborate the non-monophyly of P 
vlangalii species complex. In contrast with Guo and Wang 
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Figure 6 Bayesian skyline plots based on a mean mutation rate of 1.5x10* per site per year. The x-axis is in units of thousand years ago 
and the y-axis represents the log-scale estimated population size [(Ne*t), the product of the female effective population size and generation 
length in years]. The central line represents the median value for the population size, while the shaded area represents the 95% credibility 


limits. Lineage designations correspond to those in Figure 2. 


(2007), our results indicate that P. forsythii is not the basal 
species in the viviparous group. However, these results 
should be interpreted with caution. 

The incongruence across genes is reasonable. It is an 
inescapable biological reality that gene trees differ for a 
variety of reasons (reviewed in Maddison, 1997; Degnan 
and Rosenberg, 2009). Mitochondrial DNA sequence 
variation in animals is notoriously characterized by a 
high amount of homoplasy, i. e., genealogical conflict 
between sites, which decreased the efficiency of mtDNA 
genes for phylogenetic inference (Engstrom et al., 
2004). Historically, it is assumed that the genealogical 
relationships within a single gene mirrored those that 
actually occurred in the formation of species. However, 
it is now widely acknowledged that such direct 
correspondence between the evolutionary history of genes 
and of species divergence does not necessarily exist. In 
cases of rapid speciation, for example, it may actually 
be more common for gene histories to disagree with one 
another than to match. Given that conflicting genealogical 
histories often exist in different genes throughout the 
genome, it is imperative to take advantage of the direct 
estimation of species trees, 1. e., methods for multilocus 
species tree inference (Knowles, 2009; Knowles and 


Kubatko, 2010) to infer a robust phylogeny of the 
viviparous group of Phrynocephalus. 


4.2 Taxonomic implications As can be seen from Figure 
2, the populations around Qinghai Lake and those from 
Tianzhu and Guide together form an independently 
evolving lineage, which seem to be the sister taxon to the 
rest of viviparous species. In support of several previous 
studies (Guo and Wang, 2007; Jin, 2008; Nobel et al., 
2010; Wang et al., 2002), we recognize that P. putjatia is 
a valid species. In addition, we propose that this species 
also comprises the populations from North of Qinghai 
Lake, for example, the populations from Haiyan and 
Gangcha in this study. 

Lineage Al consisted of the populations from the 
source region of the Yellow River, including Hongyuan, 
Xiaman, Maqu and Maduo (Sites 1–7, 22 in Figure 
1). This range covers the type locality of P. vlangalii 
hongyuanensis. However, Bedriaga (1909) assigned 
populations of P. vlangalii in the Bailong Jiang region, 
to P. v. pylzowi based on morphological differences. In 
support of Jin (2008), we agree that P. v. hongyuanensis 1s 
a synonym of P. v. pylzowi. 


4.3 Development of divergent lineages within P. 


Мо. 1 


vlangalii species complex Our current estimates might 
not fully reflect the actual diversification time-scales of 
the P. vlangalii species complex because they are based 
exclusively on mitochondrial data. However, these 
estimates provide a preliminary timeframe to trace the 
intraspecific differentiation of the maternal lineages of 
P. vlangalii/P. putjataia. Two alternative hypotheses 
are suggested to explain the development of divergent 
lineages within a species: the extremely large ancestral 
population size; or allopatric differentiation in the past 
(Avise, 2000). 

Under the first scenario, we followed Luikart et al. 
(2001) to compute the possible effective size (Ne) of an 
ancestral population necessary to maintain three divergent 
mtDNA lineages for P. vlangalii. To perform this 
calculation, we assumed the coalescent time of lineages 
Al, A2 and A3 in the ancestral population to be around 
950 000 years ago (Figure 4), or 316 666 generations 
of 3 years. The expectation of this time is 2Ne (1— 1/k) 
generations, where К is 105 mtND4-tRNA""" sequences, 
implying a Ne value of 159 855. The population census 
size would have to be far larger than this estimated Ne 
(Avise, 2000). It is unlikely that a population of this size 
was present at any time. Such a large ancestral population 
is also unlikely for any other toad-headed lizards. 
Similarly, under the first scenario, we may estimate a 
Ne value of 399 975 for the ancestral population of P. 
putjatia. It is also unlikely that a population of this size 
was present at any time for P. ршјапа. 

Therefore the second scenario is more likely: 
namely that the glacial climate drove the lizards into 
a number of refugia, where the populations diverged. 
The divergence time of the three lineages of P. vlangalii 
estimated here (ca. 0.95 Ma) is roughly consistent with 
the antepenultimate glaciation events (Kunlun Glaciation) 
recorded for the QTP during the middle-late Pleistocene 
(Zheng et al., 2002), which is correlated with MIS 18-16 
(0.78—0.62 Ma). During the most extensive glaciation in 
the plateau, 1. e., Kunlun Glaciation, there were glaciers 
developing in the Bayun Har and A’nyemagen Mountains 
and four big ice caps with a diameter of about 30—50 
km. In addition, the subsequent and relatively extensive 
Zhonglianggan Glaciation occurred between 500 000 
and 400 000 years аро (Owen et al., 2006; Shi et al., 
1998). Although temperatures rose during the interglacial 
stages, development of glaciers and/or extremely low 
temperatures at high elevations (> 4 500 m) throughout 
the Pleistocene may have continued to impede the gene 
flow of P. vlangalii, eventually resulting in the genetic 
divergence of P. vlangalii in the allopatric regions. 
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In addition, it is worth noting that the A’nyemaqen 
Mountains appear to have shown significant uplift in 
the middle Pleistocene during the Kunlun-Huanghe 
movement (Cui et al., 1998), although precise dates are 
not available. Thus, in support of Jin et al. (2008) and 
Wang et al. (2009), we confirm the hypothesis that uplift 
of the A'nyemaqen Mountains and glaciations since mid- 
late Pleistocene have promoted the allopartric divergence 
of P. vlangalii in the upper reaches region of the Yellow 
River. 

The Huangshui River valley 15 situated in eastern 
Qinghai, China, and between latitude 35?56'-37?38'N 
and longitude 100?35'-103?05'E. The area belongs to 
the transition zone between Qinghai-Tibetan Plateau and 
Loess Plateau with altitudes ranging from 1650 m to 4395 
m. The divergence time of the two lineages of P. putjatia 
is estimated as about 1.78 Ma (95% HPD: 0.92-3.06 
Ma), which is highly consistent with the final phase of 
the Qingzang Movement (1.7 Ma; Li et al., 1996). During 
this phase, there took place accelerated rising at the 
Xining-Huzhu region, along Huangshui River, the first- 
class tributary of Yellow River (Lu et al., 2004). There 
was a significant readjustment of the fluvial catchment 
during 1.55—1.2 million years ago in the region of paleo- 
Huangshui, resulting in the modern pattern of Huangshui 
River. There is no phrynocephalus lizards distribution 
record in the Huangshui River valley and east of Riyue 
Shan (Zhao, 1997). Tectonic movement in this region 
was prevalent during this period and is a likely cause 
of this fragmentation. Subsequently, the glacial climate 
throughout the Pleistocene may have continued to impede 
the gene flow of P. putjatia, eventually resulting in the 
genetic divergence of P. putjatia in the allopatric regions, 
e. g., north and southeast of the Qinghai Lake. 
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Appendix 1 Tetrahedral plots. (A): all sites of ND4-tRNA‘””; (B): first codon sites; (C): second 
codon sites; (D): third codon sites; and (E): tRNAs sites. The plots were obtained using the Choice of 


Sites command from the View menu in the program SeqVis v.1.5 (Ho et al., 2006). 


H28 | Ad 


86 H26 | Р.У. pyizowi CladeA 


P vlangalii 


22 H38 H A3 
63) AY054080 
85 . | AY054081 
100 [у AY054084 


AY054085 | P theobaldi 


"v 


. v. vlangalii 


8 
I 
& 
B 


77 95 87 | AY054082 
AY054083 
100 | дуо54086 
AY054087 
i Pel HM235538 
95 wM P forsythii 
HM235555 
HM235525 
92]H30 
H31 
H1 
99 54 нив 
H17 
H16 
H7 
58 H2 B1 
H3 
99 НЕ 
i Clade B 
P. putjatia На 
H20 
100 n 
H9 
65r H12 
es | H13 
H10 B2 
100 |! H11 
H14 
99 H15 
92 AY054023 P. frontalis 


AY054060 F. przewalskii 
AY054002 F. albolineatus 
НМ235556 Р. axillaris 


10 steps 


Appendix 2 Maximum parsimony consensus tree from 1 000 bootstrap replicates of Мр4- ЕМА“ by 
using PAUP*, with gaps treated as missing data. Numbers above the branch represent percent recovery 


in bootstrap analysis (1 000 pseudoreplicates). Tree length = 423, CI = 0.5816, RI = 0.9052. 


